% estimaLoopOuter.m
% Formerly ESTIMA_LOOP_PARALLEL
% This loop has been written for Paralell Computation


%% 1.Organize Output storage (Broadcast variables)

%% loop_ and parfor 
%  parfor does not like it in structure form, has trouble slicing them 
%  first three matrices are created in Initialize.m as patfor does not
%  like indexing column 1 and column 2 in the same loop 
% loop_logPrior(:,2)=-inf*ones(priorStru.numStartVals,1);
% loop_logLikel(:,2)=-inf*ones(priorStru.numStartVals,1);
% loop_logPost(:,2)=-inf*ones(priorStru.numStartVals,1);

loop_parMode=nan(numpar, priorStru.numStartVals);
loop_hessian=zeros( length(posStru.param.est), length(posStru.param.est) , priorStru.numStartVals ); 
loop_std=nan(length(posStru.param.est), priorStru.numStartVals); 
loop_exitFlag=nan(1, priorStru.numStartVals);
loop_iterations=nan(1, priorStru.numStartVals);
loop_time=nan(1, priorStru.numStartVals);

%% 2. Assign default values to ESTIMOPT if they do not exist
if ~exist('estimopt', 'var'); estimopt = []; end;
[~,estimopt]=ch_field(estimopt,'TolFun',1e-5);
[~,estimopt]=ch_field(estimopt,'MaxFuneval',200000);
[~,estimopt]=ch_field(estimopt,'MaxIter',400);
[~,estimopt]=ch_field(estimopt,'Display','iteration');
estimopt.UseParallel='always';

%% 3  Preliminary assigments and truncations before loop begins
%     Formerly known as INDMAT, this is the matrix of transformations
%     PARVECMODE is an empty vector with the calibrated values in
parestmat=zeros(length(posStru.param.est),2);

%% 4. flags.parallel==0 only

%% 5. Additional inputs used in the estimation (estimStructure)
% estimStructure.parVec=parvec;
% estimStructure.parPosEst=posStru.param.est;
% estimStructure.transMat=partrmat(posStru.param.est, :);
% estimStructure.trainVec=trainvec;
% estimStructure.ssPosEst=ssposest;
% estimStructure.flagSilent=flags.parallel;
struInOpt.transMat=prpar.transMat(posStru.param.est,:); 

tic;

%% 6. Assign Optimization function (funcOptimization)
flags.optimization=0;
if isfield(estimopt,'funcName')==false || ...
        isempty(estimopt.funcName) == true
    disp('No Optimization');
    %     estimopt.funcName='csminwel';
    %     estimopt.funcHandle=@modeFinderCsminwel;
    %     flags.optimization=1;
else
    if strcmpi(estimopt.funcName,'FMINCON')==true;
        estimopt.funcHandle=@modeFinderFmincon;
        disp('Using fmincon for Optimization');
        flags.optimization=1;
        flags.transform=0; 
    else
        estimopt.funcName='csminwel';
        estimopt.funcHandle=@modeFinderCsminwel;
        flags.optimization=1;
        flags.transform=1; 
    end
end
if flags.optimization~=1
    error('Could Not Assign an Optimization Routine')
end

%% Ensure using the right calibrated parameters 
parStru.param(posStru.param.cal)=parStru.cal; 
%% 7. Assign Loop depending on the estimation case 
switch flags.parallel 
    case 0 
        modeFinderLoopSingle
    case 1 
        modeFinderLoopParallel 
end 

%% 9. Create Optimat, matrix of Iterations | Exit Flag and Mat_
optimat=[loop_iterations; loop_exitFlag; loop_time];
dispaj('Total time in minutes=',sum(optimat(3,:)));

%% 8. Save Workspace 


loopStru.logPost=loop_logPost;
loopStru.logLikel=loop_logLikel;
loopStru.logPrior=loop_logPrior;
loopStru.parMode=loop_parMode;
loopStru.std=loop_std;
loopStru.hessian=loop_hessian;
loopStru.time=loop_time;
loopStru.exitFlag=loop_exitFlag;
loopStru.iterations=loop_iterations;
clear ii ZZ add2Func pargrid ans two; 
clear loop_*; 

try
    cd(location.savePath);
    save workspace;
    cd(cucd);
    disp('Saved Workspace'); 
catch
    disp('Could not save workspace');
end


%% 1. Assign Flag Structure 
flags.hessianDone=0; 
flags.spectrumDone=0; 
flags.tableDone=0; 
flags.handleChanged=0; 
disp('Assigned flags'); 